library(stargazer)
library(magrittr) # for %$% pipes
library(tidyverse)
library(haven)

# Read the 'titanic.dta' data and create new variables 'd' and 's'
titanic <- read_data("titanic.dta") %>% 
  mutate(d = case_when(class == 1 ~ 1, TRUE ~ 0)) %>% # 'd' is 1 if 'class' is 1, 0 otherwise
  mutate(s = case_when(sex == 0 & age == 1 ~ 1,        # 's' is a categorical variable based on 'sex' and 'age'
                       sex == 0 & age == 0 ~ 2,
                       sex == 1 & age == 1 ~ 3,
                       sex == 1 & age == 0 ~ 4,
                       TRUE ~ 0))

# Calculate the mean of 'survived' for different groups
ey11 <- titanic %>% filter(s == 1 & d == 1) %$% mean(survived)
ey10 <- titanic %>% filter(s == 1 & d == 0) %$% mean(survived)
ey21 <- titanic %>% filter(s == 2 & d == 1) %$% mean(survived)
ey20 <- titanic %>% filter(s == 2 & d == 0) %$% mean(survived)
ey31 <- titanic %>% filter(s == 3 & d == 1) %$% mean(survived)
ey30 <- titanic %>% filter(s == 3 & d == 0) %$% mean(survived)
ey41 <- titanic %>% filter(s == 4 & d == 1) %$% mean(survived)
ey40 <- titanic %>% filter(s == 4 & d == 0) %$% mean(survived)

# Calculate the differences in means
diff1 = ey11 - ey10
diff2 = ey21 - ey20
diff3 = ey31 - ey30
diff4 = ey41 - ey40

# Calculate the number of observations where 'd' is 0
obs = nrow(titanic %>% filter(d == 0))

# Calculate weights for each group
wt1 <- titanic %>% filter(s == 1 & d == 0) %>% nrow() / obs
wt2 <- titanic %>% filter(s == 2 & d == 0) %>% nrow() / obs
wt3 <- titanic %>% filter(s == 3 & d == 0) %>% nrow() / obs
wt4 <- titanic %>% filter(s == 4 & d == 0) %>% nrow() / obs

# Calculate the weighted average treatment effect
wate = diff1*wt1 + diff2*wt2 + diff3*wt3 + diff4*wt4

# Put wate and sdo into a data frame
results <- data.frame(wate = wate, sdo = sdo)

# Use stargazer to create a table
stargazer(results, type = "text")

# Load the necessary libraries
library(tidyverse)
library(haven)
library(gridExtra)

# Define a function to read data from a .dta file hosted on GitHub
read_data <- function(df)
{
  # Construct the full URL to the file
  full_path <- paste("https://github.com/scunning1975/mixtape/raw/master/", df, sep = "")
  
  # Read the .dta file from the URL and return the resulting data frame
  df <- read_dta(full_path)
  return(df)
}

# Use the read_data function to read the "cps_mixtape.dta" file
# Bind the rows of the new data frame with the existing 'nsw_dw' data frame
# Create new variables in the data frame
nsw_dw_cpscontrol <- read_data("cps_mixtape.dta") %>% 
  bind_rows(nsw_dw) %>% 
  mutate(agesq = age^2,
         agecube = age^3,
         educsq = educ*educ,
         u74 = case_when(re74 == 0 ~ 1, TRUE ~ 0),
         u75 = case_when(re75 == 0 ~ 1, TRUE ~ 0),
         interaction1 = educ*re74,
         re74sq = re74^2,
         re75sq = re75^2,
         interaction2 = u74*hisp)

# Estimate a logistic regression model with 'treat' as the dependent variable
logit_nsw <- glm(treat ~ age + agesq + agecube + educ + educsq + 
                   marr + nodegree + black + hisp + re74 + re75 + u74 +
                   u75 + interaction1, family = binomial(link = "logit"), 
                 data = nsw_dw_cpscontrol)

# Add the predicted probabilities from the logistic regression model to the data frame
nsw_dw_cpscontrol <- nsw_dw_cpscontrol %>% 
  mutate(pscore = logit_nsw$fitted.values)

# Calculate the mean predicted probability for the control group
pscore_control <- nsw_dw_cpscontrol %>% 
  filter(treat == 0) %>% 
  pull(pscore) %>% 
  mean()

# Calculate the mean predicted probability for the treated group
pscore_treated <- nsw_dw_cpscontrol %>% 
  filter(treat == 1) %>% 
  pull(pscore) %>% 
  mean()

# Create a histogram of the predicted probabilities for the control group
hist_control <- nsw_dw_cpscontrol %>% 
  filter(treat == 0) %>% 
  ggplot() +
  geom_histogram(aes(x = pscore)) +
  ggtitle("Control Group")

# Create a histogram of the predicted probabilities for the treated group
hist_treated <- nsw_dw_cpscontrol %>% 
  filter(treat == 1) %>% 
  ggplot() +
  geom_histogram(aes(x = pscore)) +
  ggtitle("Treated Group")

# Arrange the two histograms in a grid
grid.arrange(hist_control, hist_treated, nrow = 2)